Figure 1 - This figures shows…
library(png)
library(maptools)
Checking rgeos availability: TRUE
library(raster)
library(gam)
Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14
library(mapproj)
Loading required package: maps
library(rgl)
setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")
par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)
Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.
domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)
[1] 77 8
| Region | Species | Start.Exploitation | Finish.Exploitation | Start.predomestication | Finish.predomestication | Start.Domestication | Finish.Domestication |
|---|---|---|---|---|---|---|---|
| Southwest asia | Wheat | 12.00 | 11.25 | 11.25 | 11.00 | 11.00 | 9.00 |
| Southwest asia | Barley | 12.00 | 11.25 | 11.25 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Lentil | 12.00 | 11.00 | 11.00 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Pea | 11.50 | 11.00 | 11.00 | 10.00 | 10.00 | 8.50 |
| Southwest asia | Chickpea | 11.00 | 10.50 | 10.50 | 10.25 | 10.25 | 8.25 |
| Southwest asia | Broadbean | NA | NA | NA | NA | 10.50 | NA |
| Southwest asia | Flax | 12.00 | 9.50 | NA | NA | 9.50 | NA |
| Southwest asia | Olive | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| Southwest asia | Sheep | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Goat | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Pig | 12.00 | 11.50 | 11.50 | 9.75 | 10.25 | 9.00 |
| Southwest asia | Cattle, taurine | 11.50 | 10.50 | 10.50 | 10.25 | 10.25 | 8.00 |
| Southwest asia | Cat | NA | NA | 10.50 | 4.00 | 4.00 | NA |
| South Asia | Tree cotton | 8.50 | 4.50 | NA | NA | 4.50 | NA |
| South Asia | Rice | 8.00 | 5.00 | 5.00 | 4.00 | 4.00 | 2.50 |
| South Asia | Little millet | NA | NA | NA | NA | 4.50 | NA |
| South Asia | Browntop millet | NA | NA | NA | NA | 4.00 | NA |
| South Asia | Mungbean | NA | NA | 4.50 | 3.50 | 3.50 | 3.00 |
| South Asia | Pigeonpea | NA | NA | NA | NA | 3.50 | NA |
| South Asia | Zebu cattle | 9.00 | 8.00 | NA | NA | 8.00 | 6.50 |
| South Asia | Water buffalo | 6.00 | 4.50 | NA | NA | 4.50 | NA |
| East Asia | Broomcorn Millet | 10.00 | 8.00 | NA | NA | 8.00 | NA |
| East Asia | Foxtail millet | 11.50 | 7.50 | NA | NA | 7.50 | NA |
| East Asia | Rice | 10.00 | 8.00 | 8.00 | 7.50 | 7.50 | 5.00 |
| East Asia | Soybean | 8.50 | 5.50 | NA | NA | 5.50 | 4.00 |
| East Asia | Ramie | NA | NA | NA | NA | 5.25 | NA |
| East Asia | Melon | 7.00 | 4.00 | NA | NA | 4.00 | 3.75 |
| East Asia | Pig | 12.00 | 8.50 | NA | NA | 8.50 | 6.00 |
| East Asia | Silkworm | 7.00 | 5.25 | NA | NA | 5.25 | NA |
| East Asia | Yak | NA | NA | NA | NA | 4.25 | NA |
| East Asia | Horse | 7.50 | 6.75 | 6.75 | 5.50 | 5.50 | 4.00 |
| East Asia | Bactrian Camel | NA | NA | NA | NA | 4.50 | NA |
| East Asia | Duck | 2.50 | 1.00 | NA | NA | 1.00 | NA |
| East Asia | Chicken | 6.00 | 4.00 | NA | NA | 4.00 | NA |
| New Guinea | Banana | 10.00 | 7.00 | 7.00 | 4.00 | 4.00 | NA |
| New Guinea | Taro | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| New Guinea | Yam | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| Africa and Arabia | Date palm | 7.00 | 6.00 | NA | NA | 5.00 | NA |
| Africa and Arabia | Sorghum | 8.00 | 4.00 | NA | NA | 4.00 | NA |
| Africa and Arabia | Pearl millet | NA | NA | NA | NA | 4.50 | 3.50 |
| Africa and Arabia | Fonio | NA | NA | NA | NA | 2.50 | NA |
| Africa and Arabia | Cowpea | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Hyacinth bean | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Rice | 3.50 | 2.00 | NA | NA | 2.00 | NA |
| Africa and Arabia | Oil palm | 9.25 | 3.50 | NA | NA | 3.50 | NA |
| Africa and Arabia | Cattle | NA | NA | 9.00 | 7.75 | 7.75 | 6.50 |
| Africa and Arabia | Donkey | 9.00 | 5.50 | NA | NA | 5.50 | 3.50 |
| Africa and Arabia | Dromedary camel | 6.50 | 3.00 | NA | NA | 3.00 | NA |
| Africa and Arabia | Guinea fowl | NA | NA | 2.50 | 1.50 | 1.50 | NA |
| North America | Squash | 6.50 | 5.00 | NA | NA | 5.00 | NA |
| North America | Sunflower | 6.00 | 4.75 | NA | NA | 4.00 | NA |
| North America | Sumpweed | 6.00 | 4.50 | NA | NA | 4.00 | NA |
| North America | Pitseed goosefoot | 4.75 | 3.75 | NA | NA | 3.75 | NA |
| Meso-america | Squash (pepo) | NA | NA | NA | NA | 10.00 | 9.50 |
| Meso-america | Maize | 10.00 | 9.00 | NA | NA | 9.00 | NA |
| Meso-america | Foxtail millet-grass | NA | NA | NA | NA | 6.00 | 4.00 |
| Meso-america | Common bean | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Avocado | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Chile pepper | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Turkey | NA | NA | NA | NA | 2.00 | NA |
| South America | Chili pepper | NA | NA | NA | NA | 6.00 | NA |
| South America | Peanut | NA | NA | NA | NA | 5.00 | NA |
| South America | Cotton | NA | NA | NA | NA | 6.00 | NA |
| South America | Coca | NA | NA | NA | NA | 8.00 | NA |
| South America | Now-minor root crops (arrowroot, leren) | NA | NA | NA | NA | 9.00 | NA |
| South America | Squash | NA | NA | NA | NA | 10.00 | NA |
| South America | Common bean | NA | NA | NA | NA | 5.00 | NA |
| South America | Lima bean | NA | NA | 8.25 | NA | 6.00 | NA |
| South America | Monioc | NA | NA | NA | NA | 7.00 | NA |
| South America | Sweet potato | NA | NA | NA | NA | 5.00 | NA |
| South America | White potato | 7.00 | 4.50 | NA | NA | 4.00 | NA |
| South America | Quinoa | 5.00 | NA | NA | NA | 3.50 | NA |
| South America | Yam | NA | NA | NA | NA | 5.50 | NA |
| South America | Llama | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| South America | Alpaca | 10.00 | 5.00 | NA | NA | 5.00 | NA |
| South America | Guinea pig | NA | NA | NA | NA | 5.00 | NA |
| South America | Muscovy Duck | NA | NA | NA | NA | 4.00 | NA |
par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset" )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)
hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)
par(mfrow=c(2,3), mar=c(4,4,2,0))
dim(domestication_times)
[1] 77 8
specific_dates <- domestication_times[,3:8]
for(i in c(1, 3, 5, 2, 4, 6)){
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main= names(specific_dates)[i])
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
}
I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.
h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")
I think these are best described by a cummulative distribution, showing how they accumulate over time.
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(levels(domestication_times$Region)[ type_number])
print(match)
print(j)
}
[1] "Africa and Arabia"
[1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
[1] "East Asia"
[1] 10.00 11.50 10.00 8.50 5.25 7.00 12.00 7.00 4.25 7.50 4.50 2.50 6.00
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
[1] "Meso-america"
[1] 10 10 6 3 3 3 2
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
[1] "South America"
[1] 6.0 5.0 6.0 8.0 9.0 10.0 5.0 6.0 7.0 5.0 7.0 5.0 5.5 10.0 10.0 5.0 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
[1] "Southwest asia"
[1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
break_one <- maxer
break_two <- maxer - quantile(j)[2]
polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
abline(h = 25, col="limegreen", lwd=2)
}
Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.
library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)
map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
# need to reproject
This is obviously a bad projection fit right now.
origins <- readShapePoly('Origins_updated.shp')
proj4string(origins) <- CRS("+proj=longlat +datum=WGS84")
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#proj <- CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
#origins.ea <- spTransform(origins, proj)
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 14)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA Fertile_Cresc Chinese_loess New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
20 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Fertile_Cresc ... West Africa T
origins_subset$name
NULL
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
LW <- 4
lines(x = c(-100 , -130), y = c(20 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Mesoamerica
lines(x = c(-80 , -70), y = c(0 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #NW_Lowland_SA
lines(x = c(-74 , -5 ), y = c(5 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #N_Lowland_SA
lines(x = c(40 , 30), y = c(36 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Fertile_Cresc
lines(x = c(110 , 70), y = c(40 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Chinese_loess
lines(x = c(142 , 160), y = c(-5 , 90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #New_Guinea
lines(x = c(-85 , -130), y = c(33 , -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #E_North_Ameri
lines(x = c(-68 , -75), y = c(-25 , -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #C/S_Andes
lines(x = c(-10 , -20), y = c( 15, -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #W_African_Sav
lines(x = c(25 , 40), y = c(9 , -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Sudanic_Savan
lines(x = c(87 , 100), y = c(20 , -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Ganges_E_Indi
lines(x = c(120 , 160), y = c(30 , -90), lty=2, col=adjustcolor("white", alpha=1), lwd=LW) #Lower-MiddleY
dev.off()
null device
1
#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL
class : RasterStack
dimensions : 288, 720, 207360, 18 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : fourK, fiveK, sixK, sevenK, eightK, nineK, tenK, elevenK, twelveK, thirteenK, fourteenK, fifteenK, sixteenK, seventeenK, eighteenK, ...
min values : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ...
max values : 2.051069, 2.013452, 2.142908, 1.888403, 1.863014, 1.880628, 1.650615, 1.678033, 1.697732, 1.499115, 1.517264, 1.443677, 1.464867, 1.453581, 1.436394, ...
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r
class : RasterLayer
dimensions : 288, 720, 207360 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : fourK
values : 5.611358e-07, 2.051069 (min, max)
We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.
We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.
# Read the polygons
library(rgdal)
rgdal: version: 1.2-5, (SVN revision 648)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
Path to GDAL shared files:
Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
Path to PROJ.4 shared files: (autodetected)
WARNING: no proj_defs.dat in PROJ.4 shared files
Linking to sp version: 1.2-3
getwd()
[1] "/Users/Ty/Desktop/Botero postdoc 2016/Human density and the origins of agriculture"
origins <- readShapePoly('Origins_updated.shp')
proj4string(origins)
[1] NA
# Extract data
library(raster)
e <- extent(-180, 180, -60, 84)
all_cells <- extract(r, e, cellnumber = TRUE)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
length(all_cells)
[1] 414720
for(i in 1:20){
all_cells <- all_cells[-which(per.origin[[i]][,1] %in% all_cells[,1]), ]
}
length(all_cells)
[1] 408656
names(per.origin) <- origins@data[, 1]
str(all_cells)
num [1:204328, 1:2] 3033 3034 3035 3036 3037 ...
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:2] "cell" "value"
str(per.origin)
List of 20
$ W_African_Sav: num [1:309, 1:2] 92506 92507 92508 92509 92510 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sudanic_Savan: num [1:306, 1:2] 99050 99051 99052 99053 99054 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ West Africa T: num [1:427, 1:2] 102609 102610 102611 103326 103327 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ethipian plat: num [1:275, 1:2] 106281 106282 106283 106998 106999 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Fertile_Cresc: num [1:195, 1:2] 67404 67405 67406 68119 68120 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ E_North_Ameri: num [1:170, 1:2] 64270 64271 64984 64985 64986 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ New_Guinea : num [1:24, 1:2] 127360 127361 127362 128080 128081 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Mesoamerica : num [1:73, 1:2] 93032 93033 93034 93035 93036 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ N_Lowland_SA : num [1:39, 1:2] 110373 110374 110375 111092 111093 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NW_Lowland_SA: num [1:24, 1:2] 123319 123320 123321 123322 123323 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sava_W_India : num [1:34, 1:2] 83312 84031 84032 84749 84750 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ S_India : num [1:18, 1:2] 99153 99154 99872 99873 99874 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ganges_E_Indi: num [1:92, 1:2] 85494 85495 85496 85497 85498 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Lower-MiddleY: num [1:84, 1:2] 72598 72599 73318 73319 73320 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Japanese : num [1:36, 1:2] 59681 59682 60401 61121 61843 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ W_Yuman_E_Tib: num [1:131, 1:2] 72547 72548 72549 73267 73268 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ South trop ch: num [1:178, 1:2] 84818 84819 84820 84821 84822 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Chinese_loess: num [1:258, 1:2] 62493 62494 62495 62496 62497 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Southwes amaz: num [1:194, 1:2] 137758 137759 137760 137761 137762 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ C/S_Andes : num [1:165, 1:2] 137736 137737 137738 137739 137740 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
origin_vectors <- rep(NA, 3)
for(h in 1:20){
originI <- Pop[per.origin[[h]][, 1], ]
cell_vector <- as.vector(per.origin[[h]][, 1])
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- as.vector(originI)
all_vectors_pre <- cbind(rep(names(per.origin)[h], length(x_value_vector)),x_value_vector, y_value_vector, cell_vector)
origin_vectors <- rbind(origin_vectors, all_vectors_pre)
}
#names(origin_vectors)[1] <- "location"
origin_vectors <- origin_vectors[-1,]
origin_vectors
levels(all_vectors[,1])
[1] "C/S_Andes" "Chinese_loess" "E_North_Ameri" "Ethipian plat" "Fertile_Cresc"
[6] "Ganges_E_Indi" "Japanese" "Lower-MiddleY" "Mesoamerica" "N_Lowland_SA"
[11] "New_Guinea" "not_origin" "NW_Lowland_SA" "S_India" "Sava_W_India"
[16] "South trop ch" "Southwes amaz" "Sudanic_Savan" "W_African_Sav" "W_Yuman_E_Tib"
[21] "West Africa T"
all_vectors <- cbind(all_vectors, scale(as.numeric(all_vectors[,3])))
colnames(all_vectors) <- c("location_name", "x_values", "density_values", "cell_ID", "scaled_density_values")
all_vectors <- all_vectors[, c(1,2,3,5,4)]
all_vectors